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We study non-locally coupled noisy integrate-and-fire neurons with the Fokker- 
Planck equation. A propagating pulse state and a wavy state appear as a phase 
transition from an asynchronous state. We also find a solution in which travel- 
ing pulses are emitted periodically from a pacemaker region. 
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Coherent oscillations are observed in neural systems such as the visual cortex 
and the hippocampus. The synchronization of the oscillators is considered to 
play important roles in neural information processing |T] . There are mainly two 
viewpoints in the research of the oscillatory activity in neural systems. In the 
first viewpoint, the activity of each neuron is expressed by the firing rate, and 
the coherent oscillation appears owing to the interaction of the excitatory and 
inhibitory neurons. Wilson-Cohen and Amari found first oscillatory behavior 
theoretically in interacting neurons 12 Ej. Recently, Robinson, Rennie and Rowe 
proposed a more elaborate model to explain various EEG rhythms and epileptic 
seizures If the spatial freedom is taken into consideration, the excitation 
wave can propagate. Wilson-Cowan performed numerical simulations of two 
layers of excitable neurons and inhibitory neurons [S3- In the second viewpoint, 
each neuron is regarded as an oscillator. Coherent oscillation appears as the 
global synchronization of the coupled oscillators. The global synchronization in 
general coupled oscillators was first studied by Winfree "B" . Kuramoto proposed 
a globally coupled phase oscillator model as a solvable model for the global syn- 
chronization |3. The leaky- integrate-fire model is one of the simplest models for 
a single neuron and often used to study dynamical behaviors of neural networks. 
Each neuron receives an input via synaptic connections from other neurons and 
it fires when the input goes over a threshold and sends out impulses to other 
neurons. In that sense, the coupling is instantaneous, and then the model is 
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called pulse-coupled oscillators. MiroUo and Strogatz studied a globally coupled 
system of the integrate-and-fire neurons, and showed that perfect synchroniza- 
tion occurs in a finite time |S| . The synchronization of pulse coupled oscillators 
has been studied in deterministic systems by many researchers [5] ^1 ^2 ■ If 
each oscillator's behavior is stochastic, the model is generalized to a noisy phase 
oscillator model and a noisy integrate-and-fire model. In the stochastic system, 
the coherent oscillation appears as an analogue of the phase transition in the 
statistical mechanics. Globally coupled noisy phase oscillators were studied in 
^1 El EI El J and globally coupled noisy integrate-and-fire model were studied 
in ^1 El El ■ The globally coupled system is a useful model for the detailed 
analyses, however, local or non-local interactions are more plausible, since neu- 
rons interact with other neurons via long axons or gap junctions. The non-locally 
coupled system of the deterministic integrate-and-fire neurons was also studied 
In this paper, we study a non-locally coupled noisy integrate-and-fire model 
with the direct numerical simulation of the Fokker-Planck equation. 
The equation for a noisy integrate-and-fire neuron is written as 

fjnf 

^ = i-bx+io+m, (1) 

where a; is a variable corresponding to the membrane potential, 6 is a positive 
parameter, Iq denotes an external input, and ^{t) is the Gaussian white noise 
satisfying {£,{t)£,{t')) = 2D5{t — t'). If x reaches a threshold 1, x jumps back to 
0. li b < 1 + Iq, each neuron fires spontaneously. The Fokker-Planck equation 
for the Langevin equation (f ) is 

dP d d'^P 

^ = -£(1 - + ^o)P{^) + + d{x)Mt), (2) 

where Jo{t) = —D{dP/dx)x=i is the firing rate. The stationary distribution 
Po{x) for the Fokker-Planck equation (2) is written as 1201 

Po{x) = p^(o)e{(--(i/2)fcx^}/o^ forx<0. 



Po(0)e 



{{ax-{l/2)bx'^}/D 



^{-az+{l/2)bz^}/D^^ 

^ 

^{-az+{l/2)bz^}/Dj^^ 



for < x < I, 
(3) 



where a = 1 + /q and P{0) is determined from the normalization condition 
I-oc Po{x)dx = I. The firing rate Jq is determined as Jq = DPq{Q)/ e'^-'^'^+^^/^)b^^}/D^^_ 
We have performed direct numerical simulation of Eq. (2) with the finite dif- 
ference method with Ax = 0.0002 and At = 2.5 x 10^''', and checked that the 
stationary probability distribution (3) is successfully obtained. 

We assume a non-locally coupled system composed of the noisy integrate- 
and-fire neurons. Each neuron interacts with other neurons via synaptic con- 
nections. Time delay exists generally for the synaptic connections. A model 
equation of the interacting noisy integrate-and-fire neurons is written as 

^^l-hx,+h+U{t), (4) 
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Figure 1: (a) Linear growth rate A as a function of wavenumber k of the uniform 
solution for gi{y) = 1.2(1.5 exp(—4|y|) — 0.1) at D = 0.025. The cross mark on 
the vertical line indicates the linear growth rate for A; = 0. (b) Linear growth 
rate A at fc = 27r/L as a function of D. (c) Time evolution of the profiles of the 
firing rate Joiy). 



where Xi denotes the dimensionless membrane potential for the iih neuron, 
^i(t) denotes the noise term which is assumed to be mutually independent, i.e., 
{^i{t)£,j{t')) — 2DSijS{t — t'), and li is the input to the iih neuron by the mutual 
interaction. The input li to the ith neuron from the other neurons is given by 



where tj^ is the time of the kth firing for the jth neuron, gij denotes the inter- 
action strength from the jth neuron to the jth neuron, and r denotes a decay 
constant. The sum is taken only for t > tj,. The effect of the firing of the jth 
neuron to the ith neuron decays continuously with t. If r — > 0, the coupling 
becomes instantaneous. Equation (5) is equivalent to 

^ = -{I^-J2T.9^,At-tl)}/r- (6) 

j k 

If there are infinitely many neurons at each position y, we can define the num- 
ber density of neurons with membrane potential x clearly at each position. 
The number density is expressed as n{x,y,t) at position y and time t. The 
non-locally coupled system can be studied with a mean-field approach. In the 
mean-field approach, the number density is proportional to the probability dis- 
tribution P{x, y, t) for the probability variable x. The average value of S{t — tj,) 
expresses the average firing rate at time t at the position y. It is expressed as 
jQ{y,t) = —D{dn/dx)x=i- The number density n{x,y,t) therefore obeys the 
Fokker-Planck type equation: 

^^^^ = -^(l-bx + I{y,t))n{x,y)+D^+6{x)Jo{y,t), 
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Figure 2: Peak-peak amplitude A of I{y,t) as a function of D for (a) gi{y) = 
1.2{1.5exp(-4|y|) - 0.1}, (b) = 1.5 exp(-4|y|) - 0.1 and (c) gi{y) = 

1.3{1.5exp(-4|y|) -0.1}. 

= -{Iiy,t)~Jiy,t)}/r, 
Jiy.t) = / g{y,y')Jo{y\t)dy', (7) 



where g{y,y') is the coupling strength from the neuron located at y' to the one 
at y, and /(y), Jo{y) ai'e respectively the input and the firing rate for the neuron 
at y. 

We have assumed that the time delay for the signal to transmit between y' 
and y can be neglected and g(y,y') depends only on the distance \y — y'\, i.e., 
yiUj y') = 9{\y ~ y'\)- As two simple examples of the non-local coupling, we use 
9i{y,y') = cexp{-K\y-y'\)~d and g2{y,y') = cexp(-/t|2/-2/'|) -dexp(-K'|?/- 
y'\). These forms of the coupling imply that the interaction is excitable locally, 
but the interaction strength decreases with the distance \y — y'\, and it be- 
comes inhibitory when \y — y'\ is large. This Mexican-hat type of coupling was 
used in several neural models |23 , especially to study the competitive dynamics 
in neural systems. Although two layer models of excitatory neuron layer and 
inhibitory neuron layer may be more realistic, we consider the above simpler 
one-layer model. The inhibitory interaction approaches a constant value —d for 
the coupling gi, and for the couphng 32 • The system size is assumed to be 
L = 10 as a simple example, and the periodic boundary conditions for the space 
variable y are imposed. We choose the damping constants k and k', as the 
exponential function decays to almost for the distance \y — y'\ ^ L. Therefore, 
the dynamical behaviors do not depend on the system size L qualitatively in 
the second model. But the dynamical behaviors depend on the system size L in 
the first model, because the range of the inhibitory interaction is infinite in the 
model. 

There is a stationary and uniform solution n(x, y, t) — no{x) and I{y, t) = Iq 
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Figure 3: (a) Linear growth rate A as a function oi k a.t D — 0.01 for g2{y) = 
1.2{1.5exp(— 4|y|) — 0.4exp(— (b) Peak-peak amplitudes A oi I{y,t) as a 
function of D. (c) Time evolution of the firing rate Jo{y) for the pulse train 
state at D = 0.01. 



in the non-locally coupled equation. The uniform solution satisfies 
no(x) = no(0)e^'^"-(i/')''"'>/^, for a; < 0, 



no(0)e 



{ax-(l/2)bx'^}/D 



1 - 



^{-az+(l/2)bz^}/D^^ 
J 



for < a; < 1, 
(8) 



where the parameter a is determined by the self-consistent condition 

a = 1 - gQD{dno{x)/dx)x^i, 



(9) 



where go = J g{y,y')dy' . 

To study the linear stability of the stationary and uniform solution, we 
consider small deviations 6n{x, y, t) = n{x, y, t)—no{x) and SI{y, t) — I~Iq from 
the uniform solution. The small deviations can be expressed with the Fourier 
series as Sn{x,y,t) = ^ 5n/c(a;, t) exp(iA:?/) and SI{y,t) = ^ (5/fc exp(ifc?/) under 
the periodic boundary conditions, where k — 271111/ L. The perturbations Suk 
and Slk obey coupled linear equations 



d6nk{x,t) 

m 

dShjt) 
dt 



d 



bx + lQ)5nk{x,t) + 6Ikit)no{x)} 
= -{SIkit)~g'SJokit)}/T, 



D 



dx"^ 



6ix)6Joit), 
(10) 



where SJokit) = -D{dnk/dx)^^i and g' = J g{y,y')e''^^y' -y'^dy' . For L is 
sufficiently large, g' = 2ck/(k^ + fc^) — dLSk.o for the coupling gi and g' = 



5 




Figure 4: Time evolution of the firing rate Jo{y) for the wavy state with a 
pacemaker region at D = 0.01. 

2ck/{k^ + fc^) - 2dK' /{k''^ + fc^) for the coupHng 32- The stabihty of the station- 
ary state is determined by the real part of the eigenvalues of the linear equation 
(10). But, the stationary solution no(x) is a nontrivial function of x, and it is 
not so easy to obtain the eigenvalues. Here we have evaluated the real part of 
the largest eigenvalue of the linear equation for various k by direct numerical 
simulations of Eq. (10). The dynamical behavior in the long time evolution of 
Eq. (10) is approximately expressed with the largest eigenvalue A, that is, 5nk 
and 5Ik ~ e^* for i 1. We have numerically calculated the linear growth rate 
of the norm {5nkf'dx + {dlkY}^/'^ (which grows as e'-^"'^^* for i 3> 1) every 
time-interval 0.001. Since the norm grows to infinity or decays to zero in the 
natural time evolution of the linear equation, we have renormalized the vari- 
ables every time-interval 0.001, as the norm is 1 by the rescaling cSnu — > 5nk 
and c5Ik — > 6Ik with a constant c. We have regarded the average value of 
the linear growth rate of the norm as the largest eigenvalue for Eq. (10). Fig- 
ure 1(a) displays the linear growth rate A as a function of k for the coupling 
g-^[y) = 1.2{1.5exp(-4|y|) - 0.1} at D = 0.025. The other parameters h and 
r are fixed to be 6 = 0.8, t ~ 0.01. There is discontinuity at A; = for this 
coupling. The linear growth rate at A: = takes a negative value denoted by the 
cross. The uniform state is stable for the uniform perturbation with fc = 0. The 
growth rate decreases with fc, but it is positive for k < 2. Figure 1(b) displays 
the linear growth rate as a function of D for the coupling gi at k = 2ti jL. 
The uniform state is unstable for D < 0.0291. The Hopf bifurcation occurs 
for a nonzero wavenumber. Therefore, a wavy state is expected to appear for 
D < 0.0291. We have performed direct numerical simulations for this coupling 
at D = 0.028. Figure 1(c) displays a time evolution of the profiles of the fir- 
ing rate Jo{y,t). The profile of the firing rate has a pulse structure and it is 
propagating in the right direction. Since the pulse propagates one round L with 
period T — 3.23, the velocity of the traveling pulse is L/T ~ 3.1. A regular 
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limit cycle oscillation with period T is observed at each point. The directions 
depend on the initial conditions. The traveling pulse state is an ordered state in 
the non-locally coupled system. The locally excitable interaction facilitates the 
local synchronization of the firing, but the global inhibition suppresses the com- 
plete synchronization. As a result of the frustration, a traveling pulse appears. 
The pulse state is different form the traveling pulse observed in an excitable 
system, since the uniform state is unstable in our system and the pulse state is 
spontaneously generated from the stationary asynchronous state. 

The input I{y,t) to the neuron at position y exhibits regular limit cycle 
oscillation. Figure 2(a) displays the peak-peak amplitude A, which is defined as 
the maximum value minus the minimum value of /(y, t), as a function of D. The 
oscillatory state disappears at £> = 0.02985 and the traveling pulse state changes 
into the stationary and uniform state. Inversely, as D is decreased, the traveling 
pulse state appears spontaneously from the stationary state at D = 0.0291, 
which is the critical value obtained from the linear stability analysis. That is, 
the phase transition is weakly subcritical for this coupling. We have changed the 
coupling function as gi{y) = a(1.5 exp(— 4|y|) — 0.1) with a free parameter a, 
and studied the phase transition at two other values of a = 1 and a= 1.3. The 
critical values by the linear stability analysis are Dc = 0.0232 for a = 1 and 
Dc = 0.0320 for a = 1.3. Figures 2(b) and (c) display the peak-peak amplitude 
A of I{y,t) as a function of D by the direct numerical simulation of Eq. (7). 
The bifurcation is supercritical for a = 1 and it is subcritical for a = 1.3. The 
parameter range AD = 0.02 of the hysteresis region is larger for a = 1.3 than 
the parameter range AD = 0.0075 for a = 1.2. There is a transition from the 
supercritical bifurcation to the subcritical bifurcation at a critical value slightly 
smaller than a = 1.2. 

As a second example, we consider a non-locally coupled system with the 
coupling function g2{y) = 1-8 exp(— 4|y|) — 0.48 exp(— |y|). Figure 3(a) displays 
the linear growth rate A for the stationary and uniform state as a function 
of A; at D = 0.01. The linear growth rate is a continuous function of k and 
takes a maximum at fc ~ 2. The linear growth rate takes the largest value 
at wavenumber k = 67r/10 (,i.e., wavelength L/3 ) in our finite size system 
of L = 10. The linear growth rate at k = 67r/10 takes positive values for 
D < Dc = 0.0155. We have performed direct numerical simulations for various 
D's. A wavy state with the finite wavenumber k — Gtt/IO appears in this non- 
local system for D < Dc- Figure 3(b) displays the peak-peak amplitude of I{y, t) 
as a function of D. A supercritical phase transition occurs at D 0.0155, which 
is also consistent with the linear stability analysis. Near the critical value, the 
amplitude of the oscillation is small and the wavy state seems to be sinusoidal. 
As D is decreased, the oscillation amplitude increases and the sinusoidal waves 
change into pulse trains gradually. Figure 3(c) displays the time evolution of 
the profile of the firing rate Jo{y,t) at D = 0.01. This pattern was obtained 
by decreasing D stepwise from the sinusoidal wave state near the critical point. 
The pulse number is three and it is consistent with the result of the linear 
stability analysis. The velocity of propagating pulse is w = L/{3T) = 1.67, 
where T = 2.0 is the period of oscillation at a fixed position. Figure 4 displays 
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a time evolution of a different type of wavy state. This state was obtained in 
a numerical simulation at the same parameter D = 0.01 as Fig. 3(c), starting 
from the uniform initial condition with small random perturbations. Pulses 
are created periodically near .x ^ 6 and they are propagating alternatively in 
different directions. The inversely-propagating pulses collide at a; ~ 1 and they 
disappear. Namely, there are a pacemaker region (a source region) and a sink 
region of traveling pulses in this solution. This type of wavy state including a 
pacemaker region and the simple pulse-train state are bistable. 

To summarize, we have studied the non-locally noisy integrate-and-fire model 
with the Fokker-Planck equation. Wc; have found that a traveling pulse appears 
as a result of oscillatory phase transitions. We found also a pulse-train state 
by changing the form of the interaction. The wavy states appear as a phase 
transition from a asynchronous state when the noise strength is decreased. Wc 
have investigated a one-dimensional system for the sake of simplicity of numeri- 
cal simulations, but we can generalize the model equation to a two-dimensional 
system easily. Our non-locally coupled integrate-and-fire model might be too 
simple, however, the wavy state is one of the typical dissipative structures far 
from equilibrium. Therefore, the spontaneously generated waves might be ob- 
served as some kind of brain waves also in real neural systems. 
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